{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n=============================================================\nGaussian process regression (GPR) with noise-level estimation\n=============================================================\n\nThis example illustrates that GPR with a sum-kernel including a WhiteKernel can\nestimate the noise level of data. An illustration of the\nlog-marginal-likelihood (LML) landscape shows that there exist two local\nmaxima of LML. The first corresponds to a model with a high noise level and a\nlarge length scale, which explains all variations in the data by noise. The\nsecond one has a smaller noise level and shorter length scale, which explains\nmost of the variation by the noise-free functional relationship. The second\nmodel has a higher likelihood; however, depending on the initial value for the\nhyperparameters, the gradient-based optimization might also converge to the\nhigh-noise solution. It is thus important to repeat the optimization several\ntimes for different initializations.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(__doc__)\n\n# Authors: Jan Hendrik Metzen <jhm@informatik.uni-bremen.de>\n#\n# License: BSD 3 clause\n\nimport numpy as np\n\nfrom matplotlib import pyplot as plt\nfrom matplotlib.colors import LogNorm\n\nfrom sklearn.gaussian_process import GaussianProcessRegressor\nfrom sklearn.gaussian_process.kernels import RBF, WhiteKernel\n\n\nrng = np.random.RandomState(0)\nX = rng.uniform(0, 5, 20)[:, np.newaxis]\ny = 0.5 * np.sin(3 * X[:, 0]) + rng.normal(0, 0.5, X.shape[0])\n\n# First run\nplt.figure()\nkernel = 1.0 * RBF(length_scale=100.0, length_scale_bounds=(1e-2, 1e3)) \\\n    + WhiteKernel(noise_level=1, noise_level_bounds=(1e-10, 1e+1))\ngp = GaussianProcessRegressor(kernel=kernel,\n                              alpha=0.0).fit(X, y)\nX_ = np.linspace(0, 5, 100)\ny_mean, y_cov = gp.predict(X_[:, np.newaxis], return_cov=True)\nplt.plot(X_, y_mean, 'k', lw=3, zorder=9)\nplt.fill_between(X_, y_mean - np.sqrt(np.diag(y_cov)),\n                 y_mean + np.sqrt(np.diag(y_cov)),\n                 alpha=0.5, color='k')\nplt.plot(X_, 0.5*np.sin(3*X_), 'r', lw=3, zorder=9)\nplt.scatter(X[:, 0], y, c='r', s=50, zorder=10, edgecolors=(0, 0, 0))\nplt.title(\"Initial: %s\\nOptimum: %s\\nLog-Marginal-Likelihood: %s\"\n          % (kernel, gp.kernel_,\n             gp.log_marginal_likelihood(gp.kernel_.theta)))\nplt.tight_layout()\n\n# Second run\nplt.figure()\nkernel = 1.0 * RBF(length_scale=1.0, length_scale_bounds=(1e-2, 1e3)) \\\n    + WhiteKernel(noise_level=1e-5, noise_level_bounds=(1e-10, 1e+1))\ngp = GaussianProcessRegressor(kernel=kernel,\n                              alpha=0.0).fit(X, y)\nX_ = np.linspace(0, 5, 100)\ny_mean, y_cov = gp.predict(X_[:, np.newaxis], return_cov=True)\nplt.plot(X_, y_mean, 'k', lw=3, zorder=9)\nplt.fill_between(X_, y_mean - np.sqrt(np.diag(y_cov)),\n                 y_mean + np.sqrt(np.diag(y_cov)),\n                 alpha=0.5, color='k')\nplt.plot(X_, 0.5*np.sin(3*X_), 'r', lw=3, zorder=9)\nplt.scatter(X[:, 0], y, c='r', s=50, zorder=10, edgecolors=(0, 0, 0))\nplt.title(\"Initial: %s\\nOptimum: %s\\nLog-Marginal-Likelihood: %s\"\n          % (kernel, gp.kernel_,\n             gp.log_marginal_likelihood(gp.kernel_.theta)))\nplt.tight_layout()\n\n# Plot LML landscape\nplt.figure()\ntheta0 = np.logspace(-2, 3, 49)\ntheta1 = np.logspace(-2, 0, 50)\nTheta0, Theta1 = np.meshgrid(theta0, theta1)\nLML = [[gp.log_marginal_likelihood(np.log([0.36, Theta0[i, j], Theta1[i, j]]))\n        for i in range(Theta0.shape[0])] for j in range(Theta0.shape[1])]\nLML = np.array(LML).T\n\nvmin, vmax = (-LML).min(), (-LML).max()\nvmax = 50\nlevel = np.around(np.logspace(np.log10(vmin), np.log10(vmax), 50), decimals=1)\nplt.contour(Theta0, Theta1, -LML,\n            levels=level, norm=LogNorm(vmin=vmin, vmax=vmax))\nplt.colorbar()\nplt.xscale(\"log\")\nplt.yscale(\"log\")\nplt.xlabel(\"Length-scale\")\nplt.ylabel(\"Noise-level\")\nplt.title(\"Log-marginal-likelihood\")\nplt.tight_layout()\n\nplt.show()"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.6.9"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}